Magnetostatics solutionΒΆ
[2]:
%%capture
%run current_density.ipynb
[3]:
##############################################################################
# Tree/Cotree gauging
##############################################################################
MESH = pde.mesh3.netgen(geoOCCmesh)
R = pde.tools.tree_cotree_gauge(MESH)
[4]:
R
[4]:
<27728x23597 sparse matrix of type '<class 'numpy.float64'>'
with 23597 stored elements in Compressed Sparse Column format>
[5]:
##############################################################################
# Assembly
##############################################################################
order = 1
phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
D = pde.int.assemble3(MESH, order = order)
# R0, RSS = pde.h1.assembleR3(MESH, space = 'P1', faces = 'ambient_face')
M = phi_H1 @ D @ phi_H1.T
K = dphix_H1 @ D @ dphix_H1.T +\
dphiy_H1 @ D @ dphiy_H1.T +\
dphiz_H1 @ D @ dphiz_H1.T
# Kn = RSS @ K @ RSS.T
phix_Hcurl, phiy_Hcurl, phiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = order)
curlphix_Hcurl, curlphiy_Hcurl, curlphiz_Hcurl = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = order)
M_Hcurl = phix_Hcurl @ D @ phix_Hcurl.T +\
phiy_Hcurl @ D @ phiy_Hcurl.T +\
phiz_Hcurl @ D @ phiz_Hcurl.T
K_Hcurl = curlphix_Hcurl @ D @ curlphix_Hcurl.T +\
curlphiy_Hcurl @ D @ curlphiy_Hcurl.T +\
curlphiz_Hcurl @ D @ curlphiz_Hcurl.T
C_Hcurl_H1 = phix_Hcurl @ D @ dphix_H1.T +\
phiy_Hcurl @ D @ dphiy_H1.T +\
phiz_Hcurl @ D @ dphiz_H1.T
curlphix_Hcurl_P0, curlphiy_Hcurl_P0, curlphiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'K', order = 0)
phix_Hcurl_P0, phiy_Hcurl_P0, phiz_Hcurl_P0 = pde.hcurl.assemble3(MESH, space = 'N0', matrix = 'M', order = 0)
dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)
KR = R.T@K_Hcurl@R
MR = R.T@M_Hcurl@R
r = dx_x @ D @ phix_Hcurl.T +\
dy_x @ D @ phiy_Hcurl.T +\
dz_x @ D @ phiz_Hcurl.T
cholKR = chol(KR)
x = cholKR.solve_A(R.T@r)
x = R@x
ux = curlphix_Hcurl_P0.T @ x
uy = curlphiy_Hcurl_P0.T @ x
uz = curlphiz_Hcurl_P0.T @ x
C:\Users\Radu\AppData\Local\Temp\ipykernel_848\1558448441.py:47: CholmodTypeConversionWarning:
converting matrix of class csr_matrix to CSC format
[6]:
##############################################################################
# Storing to vtk
##############################################################################
grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_H1_Scalar(grid, x, 'x')
pde.tools.vtklib.add_L2_Vector(grid,ux,uy,uz,'grad_x')
pde.tools.vtklib.writeVTK(grid, 'magnetostatics_solution.vtu')
[7]:
import pyvista as pv
mesh = pv.read('magnetostatics_solution.vtu')
mesh2 = pv.read('current_density.vtu')
mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,1])
p = pv.Plotter()
threshed.set_active_scalars("x")
p.add_mesh(threshed, style='surface', color="w", opacity=0.2, label=None)
threshed.set_active_vectors("grad_x")
arrows = mesh.glyph(scale="grad_x", orient=True, tolerance=0.03, factor=1000)
p.add_mesh(arrows, color="orange")
arrows2 = mesh2.glyph(scale="grad_J", orient=True, tolerance=0.03, factor=9500.0)
p.add_mesh(arrows2, color="black")
p.camera_position = [(0, -400, 400),(0, 0, 0),(0, 0, 0)]
p.show(jupyter_backend="html")